Invasion threshold in heterogeneous metapopulation networks 
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We study the dynamics of epidemic and reaction-diffusion processes in metapopulation models 
with heterogeneous connectivity pattern. In SIR-like processes, along with the standard local epi- 
demic threshold, the system exhibits a global invasion threshold. We provide an explicit expression 
of the threshold that sets a critical value of the diffusion/mobility rate below which the epidemic 
is not able to spread to a macroscopic fraction of subpopulations. The invasion threshold is found 
to be affected by the topological fluctuations of the metapopulation network. The presented results 
provide a general framework for the understanding of the effect of travel restrictions in epidemic 
containment. 
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The role of heterogeneity has been acknowledged as 
a central question in the study of population biology of 
infectious diseases [H, [3, |1] and revamped recently with 
the evidence that a large number of real world networks 
exhibit complex topological properties 0, S 0|- These 
features, often mathematically encoded in a heavy-tailed 
probability distribution P{k) that any given node has de- 
gree k, were shown to affect the system evolution alter- 
ing the threshold behavior and the associated dynamical 
phase transition 0, H, @] ■ These studies have mainly fo- 
cused on networked systems where each node corresponds 
to a single individual and only recently the study of the 
impact of heterogeneous topologies on bosonic systems, 
where nodes can be occupied by any number of parti- 
cles, has been initiated [lO|. Examples are provided by 
reaction-diffusion systems used to model a wide range of 
phenomena in chemistry and physics 11 ll. and metapop- 
ulation epidemic models [1, [1, Il3l. Il4 , 3, l3| where 
particles represent people moving across different sub- 
populations (nodes) such as city or urban areas, and the 
reaction processes account for the local infection dynam- 



Here we analyze epidemic metapopulation models 
characterized by an infection dynamics within each node 
(or subpopulation) that follows a Susceptible-Infected- 
Removed (SIR) model. The mobility ratep of individuals 
defines the coupling process among the subpopulations. 
In the real world, the networks representing the mobility 
pattern of individuals among different subpopulations are 
in many cases highly heterogeneous 17, 18, Til [20l. [2H. 
For this reason, the connectivity pattern of the metapop- 
ulation network is described as a random graph with ar- 
bitrary degree distribution P{k). By using a mechanis- 
tic approach it is possible to show that along with the 
usual epidemic threshold condition i?o > 1 on the ba- 
sic reproductive number, the system exhibits a global 
invasion threshold setting the condition for the infec- 
tion of a macroscopic fraction of the metapopulation sys- 
tem ■ The threshold condition on Rq ensures 



the local outbreak at the subpopulation level [H, \ld^ . 
whereas the explicit expression obtained for global in- 
vasion threshold i?* > 1 provides a critical value for the 
diffusion rate p below which the epidemic cannot prop- 
agate to a relevant fraction of subpopulations. We find 
that the global invasion threshold is affected by the topo- 
logical fluctuations of the underlying network. The larger 
the network heterogeneity and the smaller is the value 
of the critical diffusion rate above which the epidemic 
may globally invade the metapopulation system. The 
present results can be generalized to more realistic diffu- 
sion and mobility schemes and provide a framework for 
the analysis of realistic metapopulation epidemic mod- 
els l23,[2i[2l[2i,[2i. 



A simplified mechanistic (i.e. microscopic in the 
epidemic terminology) approach to the metapopulation 
spreading of infectious diseases uses a markovian assump- 
tion in which at each time step the movement of indi- 
viduals is given in terms of a matrix dij that expresses 
the probability that an individual in the subpopulation 
i is traveling to the subpopulation j. Several modeling 
approaches to the large scale spreading of infectious dis- 

m 



eases m, 23, 24, 25 



use this mobility process based 
on transportation networks combined with the local evo- 
lution of the disease. The markovian character lies in 
the assumption that at each time step the same traveling 
probability applies to all individuals in the subpopulation 
without having memory of their origin. This mobility 
scheme coupled with an infection dynamics at the local 
level can be generally viewed as equivalent to classic reac- 
tion diffusion processes with no constraint on the occupa- 
tion numbers iV^ of each subpopulation. The total pop- 
ulation of the metapopulation system is A'^ = Ni and 
each individual diffuses along the edges with a diffusion 
coefficient dij that depends on the node degree, subpopu- 
lation size and/or the mobility matrix. The metapopula- 
tion system is therefore composed of a network substrate 
connecting nodes - each corresponding to a subpopula- 
tion - over which individuals diffuse. We consider that 
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each node i is connected to other ki nodes according to 
its degree resulting in a network with degree distribution 
P{k) and distribution moments (fc") — k°'P{k). 

In the following, as a simplified diffusion process we 
assume that the mobility is equivalent to a diffusion rate 
along any given link of a node with degree k simply equal 
to dkk' = p/k. This is obviously not the case in a wide 
range of real systems where the extreme heterogeneity 
of traffic is well documented and more realistic processes 
will be considered elsewhere. This simple process how- 
ever automatically generates a stationary distribution of 
occupation numbers that is better described by grouping 
subpopulations according to their degree k 
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where N is the average subpopulation size. 

In each subpopulation j the disease follows an SIR 
model and the total number of individuals is partitioned 
in the compartments Sj{t), Ij{t) and Rj{t), denoting the 
number of susceptible, infected and removed individuals 
at time t, respectively. The infection dynamics proceeds 
as follows. Each susceptible individual has a transition 
rate to the infected state expressed as f3Ij/Nj, where 
f3 is the disease transmissibility rate and Ij/Nj is the 
force of infection in the homogeneous mixing assumption. 
Analogously, each infected individual enters the removed 
compartment according to the recovery rate /x. The ba- 
sic SIR rules thus define a reaction scheme of the type 
S + I ^ 21 and / — > i?, that conserves the number of 
individuals. The SIR epidemic model is characterized by 
the reproductive number Rq = that defines the av- 
erage number of infectious individuals generated by one 
infected individual in a fully susceptible population. The 
epidemic is able to generate a number of infected indi- 
viduals larger than those who recover only if -Rq > 1; 
yielding the classic result for the epidemic threshold [if. 
If the spreading rate is not large enough to allow a re- 
productive number larger than one (i.e. (3 > fj,), the 
epidemic outbreak will quickly die out. This result is 
valid at the level of each subpopulation and holds also 
at the metapopulation level where i?o > 1 is a necessary 
condition to have the growth of the epidemic [l^ . 

The intuitive result on the subpopulation epidemic 
threshold however does not take into account the effects 
due to the finite size of subpopulations, the discrete na- 
ture of individuals and the stochastic nature of the re- 
action and diffusion processes. These effects have been 
shown to have a crucial role in the problem of resurgent 
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epidemics, extinction and eradication [2j 
Also in the present framework indeed each subpopula- 
tion may or may not transmit the infection to a neigh- 
boring subpopulation upon the condition that at least 
one infected individual is moving onto the non-infected 
subpopulations during the epidemic outbreak. Given an 



SIR model with i?o > 1, the total number of infected in- 
dividuals generated within a subpopulation and the mo- 
bility rate must be large enough to ensure the seeding 
of other subpopulations before the end of the local out- 
break [11,123]. 

As a simple example of this effect let us consider a 
metapopulation system in which the initial condition is 
provided by a single infection in a subpopulation with 
degree k and Nk individuals, given Ro > 1. In the case 
of a macroscopic outbreak in a closed population the to- 
tal number of infected individuals during the outbreak 
evolution will be equal to aN^ where a depends on the 
specific disease model and parameter values used. Each 
infected individual stays in the infectious state for an av- 
erage time equal to the inverse of the recovery rate, 
during which it can travel to the neighboring subpop- 
ulation of degree k' with rate dkk'- We can therefore 
consider that on average the number of new seeds that 
may appear into a connected subpopulation of degree k' 
during the duration of the local outbreak is given by 
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In this perspective we can provide a characterization of 
the invasion dynamics at the level of the subpopulations, 
translating epidemiological and demographic parameters 
into Levins-type metapopulation parameters of extinc- 
tion and invasion rate. Let us define D° as the num- 
ber of diseased subpopulation of degree k at generation 
0, i.e. those which are experiencing an outbreak at the 
beginning of the process. Each infected subpopulation 
will seed - during the course of the outbreak - the infec- 
tion in neighboring subpopulations defining the set i?^ 
of infected subpopulations at generation 1, and so on. 
This corresponds to a basic branching process H, H, 33j 
where the n— th generation of infected subpopulations of 
degree k is denoted D^. 

In order to describe the early stage of the subpopula- 
tions invasion dynamics we assume that the number of 
subpopulations affected by an outbreak (with Rq > 1) 
is small and we can therefore study the evolution of the 
number of diseased subpopulations by using a branching 
process approximation relating Z?^ with D^~^ . Let us 
consider a metapopulation network with degree distribu- 
tion P{k) and V subpopulations and write the number 
of subpopulations of degree k invaded at the generation 
n as: 
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This equation assumes that each infected subpopulation 
of degree k' of the {n — 1)— th generation, i^^r^, will 
seed the infection in a number [k' — 1) of subpopulations 
corresponding to the number of neighboring subpopula- 
tions k' minus the one which originally transmitted the 
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infection, the probability P{k\k') that each of the k' — 1 
not yet infected neighboring subpopulations has degree fc, 
and the probabihty to observe an outbreak in the seeded 
subpopulation i.e. (l — iJ^j''^ ''''■) [s^. The last factor stems 
from the probability of extinction Pg^t = I/-R0 of an epi- 
demic seeded with a single infectious individual [s^. In 
order to obtain an explicit result we will consider in the 
following that Rq — 1 ^ 1, thus assuming that the sys- 
tem is found to be very close to the epidemic threshold. 
In this limit we can approximate the outbreak proba- 
bility as 1 — Rq''^''''' ~ ^k'k{Ro ~ !)• The case of ho- 
mogeneous diffusion of individuals d^. = p/k with the 
stationary solution of eq. ^ for the subpopulation size 
yields Xk'k = pNaiJ,~^ / {k). In addition, we assume that 
at the early stage of the epidemic D]^~^/Vk ^ 1, and we 
consider the case of uncorrelated networks in which the 
conditional probability does not depend on the originat- 
ing node, i.e. P{k\k') = kP{k)/{k) [H, obtaining 



Dl = {Ro-l 
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By defining 0" — J^k' ^k' i^' ^ 1) last expression can 
be conveniently written in the iterative form 
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that allows the increasing of infected subpopulations and 
a global epidemic in the metapopulation process only if 



(i?o - 1) 
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> 1, 
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defining the global invasion threshold of the metapopula- 
tion system. In other words, i?* is the analogous of the 
basic reproductive number at the subpopulations level 
and is a crucial indicator in assessing the behavior of 
epidemics in metapopulation models. Its expression in- 
deed contains the probability of generating an outbreak 
in a neighbor subpopulation by means of mobility pro- 
cesses, (i?o — ^)pNa/{fj.{k)) for i?o — 1 ^ 1, times the 
factor {k"^) I {k) — 1 which also appears in the threshold 
conditions characterizing phase transitions on complex 
networks 0, IH, [3, HI. The explicit form of Eq. ® 
can be used to find the minimum mobility rate ensuring 
that on average each subpopulation can seed more than 
one neighboring subpopulations. The constant a is larger 
than zero for any i?o > 1, and in the SIR case for Rq close 
to 1 it can be approximated by a ~ 2(_Ro — l)/^o [3^ - 
yielding for the SIR model a critical mobility value pc 
below which the epidemics cannot invade the metapopu- 
lation system given by the equation 
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The above condition readily tells us that the closer to the 
epidemic threshold is the single subpopulation outbreak 




FIG. 1: Phase diagram of the metapopulation system. The 
final epidemic size is shown as a function of the local threshold 
i?o and of the diffusion probability p. 



and the larger it has to be the mobility rate in order to 
sustain the global spread into the metapopulation model. 
It is important to stress that when i?o increases, the small 
(i?o — 1) expansions are no longer valid and the invasion 
threshold is obtained only in the form of a complicate 
implicit expression. 

In addition, eq. ([7]) shows the dependance of the crit- 
ical mobility on the topological fluctuations of the mo- 
bility network. The ratio (fc)^/(fc(fc — 1)) is extremely 
small in heavy-tailed networks and it is vanishing in the 
limit of infinite network size. This implies that the het- 
erogeneity of the metapopulation network is favoring the 
global spread of epidemics by lowering the global spread- 
ing threshold. In other words, the topological fluctua- 
tions favor the subpopulation invasion and suppress the 
phase transition in the infinite size limit. This finding 
provides a theoretical framework and rationale for the 
evidence concerning the inefhcacy of travel restrictions 
in the containment of global epidemics 25, 3^. The sim- 
ple plug in of the actual numbers for modern transporta- 
tion networks, the population sizes and realistic disease 
parameters in the expression ([7]) indicates that a reduc- 
tion of one order of magnitude of the mobility is not 
enough to bring the system below the invasion thresh- 
old. While more complicate mobility schemes should be 
considered for a precise calculation, this result is setting 
the framework for the understanding of mobility effects 
in the spreading and containment of infectious diseases. 

In order to support the previous analytical finding we 
have performed an extensive set of Monte Carlo numeri- 
cal simulations of the metapopulation system. The sub- 
strate network is given by an uncorrelated complex net- 
work with P(A;) ~ gen erated with the uncorre- 
lated configuration model 37[ to avoid inherent struc- 
tural correlations. Network sizes of = 10^ and 10^ 
nodes have been considered. The dynamics proceeds in 
parallel and considers discrete time steps representing the 
unitary time scale r of the process. The reaction and dif- 
fusion rates are therefore converted into probabilities and 
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at each time step in each subpopulation j a susceptible 
individual is turned into an infectious with probability 
1 — (1 — -^-tY^ and each infectious individual is sub- 
ject to the recovery process and becomes recovered with 
probability /ir. The mobility is modeled assuming a dif- 
fusion probability for each individual along each link of 
the subpopulation of the form c?j.j.' — p/k. 
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FIG. 2: Effect of the network heterogeneity on the global 
epidemic threshold. The final fraction of diseased subpopu- 
lations D{oq)/V at the end of the global epidemic is shown 
as a function of the mobility rate p in a homogeneous and a 
heterogeneous network. 

A complete analysis of the system phase diagram is 
obtained by analyzing the behavior of the global attack 
rate R{oo) /N, defined as the total fraction of cases in the 
metapopulation system at the end of the epidemic, as a 
function of both Rq and p. Figurc[l]rcports the global at- 
tack rate surface in the p-Ro space and clearly shows the 
effect of different couplings as expressed by the value of p 
in reducing the final size of the epidemic at a given fixed 
value of Rq. The smaller the value of Rq, the higher the 
coupling needs to be in order for the virus to successfully 
invade a finite fraction of the subpopulations, in agree- 
ment with the analytic result of eq. ([7]) . This provides a 
clear illustration of the varying global invasion threshold 
as a function of the reproductive rate Rq. Furthermore, it 
is possible to study the effect of the heterogeneity of the 
metapopulation structure on the global epidemic thresh- 
old. Figure [2] shows the results obtained by comparing 
two random metapopulations networks, one with pois- 
sonian degree distribution (homogeneous network) and 
one with heavy-tailed (P{k) ~ degree distribution 

(heterogeneous neywork). Despite the two models have 
the same average degree, disease parameters, the fluctu- 
ations of the power-law network increase the value of i?, 
thus lowering the critical value of the mobility. 

The present analysis provides insights in setting a 
framework for the analysis of large scale spreading of epi- 
demics in realistic mobility networks. Furthermore, these 
results open the path to future work aimed at analyzing 
refined metapopulation infection models. 
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